1. /*
  2. Fast Fourier/Cosine/Sine Transform by T. Ooura
  3. This source file is created by extracting the neccesary part from a file
  4. "fft4g.c" packed in "fft.lzh".
  5. --------------------------------------------------------------
  6. Copyright:
  7. Copyright(C) 1996-1999 Takuya OOURA
  8. email: ooura@mmm.t.u-tokyo.ac.jp
  9. download: http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html
  10. --------------------------------------------------------------
  11. dimension :one
  12. data length :power of 2
  13. decimation :frequency
  14. radix :4, 2
  15. data :inplace
  16. table :use
  17. -------- Real DFT / Inverse of Real DFT --------
  18. [definition]
  19. <case1> RDFT
  20. R[k] = sum_j=0^n-1 a[j]*cos(2*pi*j*k/n), 0<=k<=n/2
  21. I[k] = sum_j=0^n-1 a[j]*sin(2*pi*j*k/n), 0<k<n/2
  22. <case2> IRDFT (excluding scale)
  23. a[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
  24. sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
  25. sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
  26. [usage]
  27. <case1>
  28. ip[0] = 0; // first time only
  29. rdft(n, 1, a, ip, w);
  30. <case2>
  31. ip[0] = 0; // first time only
  32. rdft(n, -1, a, ip, w);
  33. [parameters]
  34. n :data length (int)
  35. n >= 2, n = power of 2
  36. a[0...n-1] :input/output data (double *)
  37. <case1>
  38. output data
  39. a[2*k] = R[k], 0<=k<n/2
  40. a[2*k+1] = I[k], 0<k<n/2
  41. a[1] = R[n/2]
  42. <case2>
  43. input data
  44. a[2*j] = R[j], 0<=j<n/2
  45. a[2*j+1] = I[j], 0<j<n/2
  46. a[1] = R[n/2]
  47. ip[0...*] :work area for bit reversal (int *)
  48. length of ip >= 2+sqrt(n/2)
  49. strictly,
  50. length of ip >=
  51. 2+(1<<(int)(log(n/2+0.5)/log(2))/2).
  52. ip[0],ip[1] are pointers of the cos/sin table.
  53. w[0...n/2-1] :cos/sin table (double *)
  54. w[],ip[] are initialized if ip[0] == 0.
  55. [remark]
  56. Inverse of
  57. rdft(n, 1, a, ip, w);
  58. is
  59. rdft(n, -1, a, ip, w);
  60. for (j = 0; j <= n - 1; j++) {
  61. a[j] *= 2.0 / n;
  62. }
  63. .
  64. */
  65. #ifndef SN_H
  66. #include "sn.h"
  67. #endif
  68. void rdft(int n, int isgn, fftType *a, int *ip, fftType *w);
  69. void makect(int nc, int *ip, fftType *c);
  70. void makewt(int nw, int *ip, fftType *w);
  71. void bitrv2(int n, int *ip, fftType *a);
  72. void cftfsub(int n, fftType *a, fftType *w);
  73. void cftbsub(int n, fftType *a, fftType *w);
  74. void cftmdl(int n, int l, fftType *a, fftType *w);
  75. void cft1st(int n, fftType *a, fftType *w);
  76. void rftfsub(int n, fftType *a, int nc, fftType *c);
  77. void rftbsub(int n, fftType *a, int nc, fftType *c);
  78. void rdft(int n, int isgn, fftType *a, int *ip, fftType *w)
  79. {
  80. int nw, nc;
  81. fftType xi;
  82. nw = ip[0];
  83. if (n > (nw << 2)) {
  84. nw = n >> 2;
  85. makewt(nw, ip, w);
  86. }
  87. nc = ip[1];
  88. if (n > (nc << 2)) {
  89. nc = n >> 2;
  90. makect(nc, ip, w + nw);
  91. }
  92. if (isgn >= 0) {
  93. if (n > 4) {
  94. bitrv2(n, ip + 2, a);
  95. cftfsub(n, a, w);
  96. rftfsub(n, a, nc, w + nw);
  97. } else if (n == 4) {
  98. cftfsub(n, a, w);
  99. }
  100. xi = a[0] - a[1];
  101. a[0] += a[1];
  102. a[1] = xi;
  103. } else {
  104. a[1] = 0.5 * (a[0] - a[1]);
  105. a[0] -= a[1];
  106. if (n > 4) {
  107. rftbsub(n, a, nc, w + nw);
  108. bitrv2(n, ip + 2, a);
  109. cftbsub(n, a, w);
  110. } else if (n == 4) {
  111. cftfsub(n, a, w);
  112. }
  113. }
  114. }
  115. void bitrv2(int n, int *ip, fftType *a)
  116. {
  117. int j, j1, k, k1, l, m, m2;
  118. fftType xr, xi, yr, yi;
  119. ip[0] = 0;
  120. l = n;
  121. m = 1;
  122. while ((m << 3) < l) {
  123. l >>= 1;
  124. for (j = 0; j < m; j++) {
  125. ip[m + j] = ip[j] + l;
  126. }
  127. m <<= 1;
  128. }
  129. m2 = 2 * m;
  130. if ((m << 3) == l) {
  131. for (k = 0; k < m; k++) {
  132. for (j = 0; j < k; j++) {
  133. j1 = 2 * j + ip[k];
  134. k1 = 2 * k + ip[j];
  135. xr = a[j1];
  136. xi = a[j1 + 1];
  137. yr = a[k1];
  138. yi = a[k1 + 1];
  139. a[j1] = yr;
  140. a[j1 + 1] = yi;
  141. a[k1] = xr;
  142. a[k1 + 1] = xi;
  143. j1 += m2;
  144. k1 += 2 * m2;
  145. xr = a[j1];
  146. xi = a[j1 + 1];
  147. yr = a[k1];
  148. yi = a[k1 + 1];
  149. a[j1] = yr;
  150. a[j1 + 1] = yi;
  151. a[k1] = xr;
  152. a[k1 + 1] = xi;
  153. j1 += m2;
  154. k1 -= m2;
  155. xr = a[j1];
  156. xi = a[j1 + 1];
  157. yr = a[k1];
  158. yi = a[k1 + 1];
  159. a[j1] = yr;
  160. a[j1 + 1] = yi;
  161. a[k1] = xr;
  162. a[k1 + 1] = xi;
  163. j1 += m2;
  164. k1 += 2 * m2;
  165. xr = a[j1];
  166. xi = a[j1 + 1];
  167. yr = a[k1];
  168. yi = a[k1 + 1];
  169. a[j1] = yr;
  170. a[j1 + 1] = yi;
  171. a[k1] = xr;
  172. a[k1 + 1] = xi;
  173. }
  174. j1 = 2 * k + m2 + ip[k];
  175. k1 = j1 + m2;
  176. xr = a[j1];
  177. xi = a[j1 + 1];
  178. yr = a[k1];
  179. yi = a[k1 + 1];
  180. a[j1] = yr;
  181. a[j1 + 1] = yi;
  182. a[k1] = xr;
  183. a[k1 + 1] = xi;
  184. }
  185. } else {
  186. for (k = 1; k < m; k++) {
  187. for (j = 0; j < k; j++) {
  188. j1 = 2 * j + ip[k];
  189. k1 = 2 * k + ip[j];
  190. xr = a[j1];
  191. xi = a[j1 + 1];
  192. yr = a[k1];
  193. yi = a[k1 + 1];
  194. a[j1] = yr;
  195. a[j1 + 1] = yi;
  196. a[k1] = xr;
  197. a[k1 + 1] = xi;
  198. j1 += m2;
  199. k1 += m2;
  200. xr = a[j1];
  201. xi = a[j1 + 1];
  202. yr = a[k1];
  203. yi = a[k1 + 1];
  204. a[j1] = yr;
  205. a[j1 + 1] = yi;
  206. a[k1] = xr;
  207. a[k1 + 1] = xi;
  208. }
  209. }
  210. }
  211. }
  212. #if USES_SIN_COS_TABLE == 0
  213. /// Ooura's version
  214. void makewt(int nw, int *ip, fftType *w)
  215. {
  216. int j, nwh;
  217. fftType delta, x, y;
  218. ip[0] = nw;
  219. ip[1] = 1;
  220. if (nw > 2) {
  221. nwh = nw >> 1; // nwh =1024, 2048, ... 2^k
  222. // delta = atan(1.0) / nwh;
  223. delta = fftM_PI_4 / nwh;
  224. w[0] = 1;
  225. w[1] = 0;
  226. w[nwh] = fftCos(delta * nwh);
  227. w[nwh + 1] = w[nwh];
  228. if (nwh > 2) {
  229. for (j = 2; j < nwh; j += 2) {
  230. x = fftCos(delta * j);
  231. y = fftSin(delta * j);
  232. w[j] = x;
  233. w[j + 1] = y;
  234. w[nw - j] = y;
  235. w[nw - j + 1] = x;
  236. }
  237. bitrv2(nw, ip + 2, w);
  238. }
  239. }
  240. }
  241. void makect(int nc, int *ip, fftType *c)
  242. {
  243. int j, nch;
  244. fftType delta;
  245. ip[1] = nc;
  246. if (nc > 1) {
  247. nch = nc >> 1; // nch = 1024, 2048,...
  248. delta = fftM_PI_4/nch; // atan(1.0) / nch;
  249. c[0] = fftCos(delta * nch);
  250. c[nch] = 0.5 * c[0];
  251. for (j = 1; j < nch; j++) {
  252. c[j] = 0.5 * fftCos(delta * j);
  253. c[nc - j] = 0.5 * fftSin(delta * j);
  254. }
  255. }
  256. }
  257. #else
  258. void makewt(int nw, int *ip, fftType *w)
  259. {
  260. // int j, nwh;
  261. long j, nwh;
  262. // fftType delta, x, y;
  263. fftType x, y;
  264. ip[0] = nw;
  265. ip[1] = 1;
  266. if (nw > 2) {
  267. nwh = nw >> 1;
  268. // delta = atan(1.0) / nwh;
  269. // delta = fftM_PI_4 / nwh;
  270. w[0] = 1;
  271. w[1] = 0;
  272. w[nwh] = fftCosR(nwh, nwh);
  273. // w[nwh] = fftCos(delta * nwh);
  274. w[nwh + 1] = w[nwh];
  275. if (nwh > 2) {
  276. for (j = 2; j < nwh; j += 2) {
  277. x = fftCosR(j, nwh);
  278. y = fftSinR(j, nwh);
  279. // x = fftCos(delta * j);
  280. // y = fftSin(delta * j);
  281. w[j] = x;
  282. w[j + 1] = y;
  283. w[nw - j] = y;
  284. w[nw - j + 1] = x;
  285. }
  286. bitrv2(nw, ip + 2, w);
  287. }
  288. }
  289. }
  290. void makect(int nc, int *ip, fftType *c)
  291. {
  292. // int j, nch;
  293. long j, nch;
  294. // fftType delta;
  295. ip[1] = nc;
  296. if (nc > 1) {
  297. nch = nc >> 1;
  298. // delta = fftM_PI_4/nch; // atan(1.0) / nch;
  299. // c[0] = fftCos(delta * nch);
  300. c[0] = fftCosR(nch, nch);
  301. c[nch] = 0.5 * c[0];
  302. for (j = 1; j < nch; j++) {
  303. // c[j] = 0.5 * fftCos(delta * j);
  304. // c[nc - j] = 0.5 * fftSin(delta * j);
  305. c[j] = 0.5 * fftCosR(j, nch);
  306. c[nc - j] = 0.5 * fftSinR(j, nch);
  307. }
  308. }
  309. }
  310. #endif // USES_SIN_COS_TABLE
  311. void cftfsub(int n, fftType *a, fftType *w)
  312. {
  313. int j, j1, j2, j3, l;
  314. fftType x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  315. l = 2;
  316. if (n > 8) {
  317. cft1st(n, a, w);
  318. l = 8;
  319. while ((l << 2) < n) {
  320. cftmdl(n, l, a, w);
  321. l <<= 2;
  322. }
  323. }
  324. if ((l << 2) == n) {
  325. for (j = 0; j < l; j += 2) {
  326. j1 = j + l;
  327. j2 = j1 + l;
  328. j3 = j2 + l;
  329. x0r = a[j] + a[j1];
  330. x0i = a[j + 1] + a[j1 + 1];
  331. x1r = a[j] - a[j1];
  332. x1i = a[j + 1] - a[j1 + 1];
  333. x2r = a[j2] + a[j3];
  334. x2i = a[j2 + 1] + a[j3 + 1];
  335. x3r = a[j2] - a[j3];
  336. x3i = a[j2 + 1] - a[j3 + 1];
  337. a[j] = x0r + x2r;
  338. a[j + 1] = x0i + x2i;
  339. a[j2] = x0r - x2r;
  340. a[j2 + 1] = x0i - x2i;
  341. a[j1] = x1r - x3i;
  342. a[j1 + 1] = x1i + x3r;
  343. a[j3] = x1r + x3i;
  344. a[j3 + 1] = x1i - x3r;
  345. }
  346. } else {
  347. for (j = 0; j < l; j += 2) {
  348. j1 = j + l;
  349. x0r = a[j] - a[j1];
  350. x0i = a[j + 1] - a[j1 + 1];
  351. a[j] += a[j1];
  352. a[j + 1] += a[j1 + 1];
  353. a[j1] = x0r;
  354. a[j1 + 1] = x0i;
  355. }
  356. }
  357. }
  358. void cftbsub(int n, fftType *a, fftType *w)
  359. {
  360. int j, j1, j2, j3, l;
  361. fftType x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  362. l = 2;
  363. if (n > 8) {
  364. cft1st(n, a, w);
  365. l = 8;
  366. while ((l << 2) < n) {
  367. cftmdl(n, l, a, w);
  368. l <<= 2;
  369. }
  370. }
  371. if ((l << 2) == n) {
  372. for (j = 0; j < l; j += 2) {
  373. j1 = j + l;
  374. j2 = j1 + l;
  375. j3 = j2 + l;
  376. x0r = a[j] + a[j1];
  377. x0i = -a[j + 1] - a[j1 + 1];
  378. x1r = a[j] - a[j1];
  379. x1i = -a[j + 1] + a[j1 + 1];
  380. x2r = a[j2] + a[j3];
  381. x2i = a[j2 + 1] + a[j3 + 1];
  382. x3r = a[j2] - a[j3];
  383. x3i = a[j2 + 1] - a[j3 + 1];
  384. a[j] = x0r + x2r;
  385. a[j + 1] = x0i - x2i;
  386. a[j2] = x0r - x2r;
  387. a[j2 + 1] = x0i + x2i;
  388. a[j1] = x1r - x3i;
  389. a[j1 + 1] = x1i - x3r;
  390. a[j3] = x1r + x3i;
  391. a[j3 + 1] = x1i + x3r;
  392. }
  393. } else {
  394. for (j = 0; j < l; j += 2) {
  395. j1 = j + l;
  396. x0r = a[j] - a[j1];
  397. x0i = -a[j + 1] + a[j1 + 1];
  398. a[j] += a[j1];
  399. a[j + 1] = -a[j + 1] - a[j1 + 1];
  400. a[j1] = x0r;
  401. a[j1 + 1] = x0i;
  402. }
  403. }
  404. }
  405. void cft1st(int n, fftType *a, fftType *w)
  406. {
  407. int j, k1, k2;
  408. fftType wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
  409. fftType x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  410. x0r = a[0] + a[2];
  411. x0i = a[1] + a[3];
  412. x1r = a[0] - a[2];
  413. x1i = a[1] - a[3];
  414. x2r = a[4] + a[6];
  415. x2i = a[5] + a[7];
  416. x3r = a[4] - a[6];
  417. x3i = a[5] - a[7];
  418. a[0] = x0r + x2r;
  419. a[1] = x0i + x2i;
  420. a[4] = x0r - x2r;
  421. a[5] = x0i - x2i;
  422. a[2] = x1r - x3i;
  423. a[3] = x1i + x3r;
  424. a[6] = x1r + x3i;
  425. a[7] = x1i - x3r;
  426. wk1r = w[2];
  427. x0r = a[8] + a[10];
  428. x0i = a[9] + a[11];
  429. x1r = a[8] - a[10];
  430. x1i = a[9] - a[11];
  431. x2r = a[12] + a[14];
  432. x2i = a[13] + a[15];
  433. x3r = a[12] - a[14];
  434. x3i = a[13] - a[15];
  435. a[8] = x0r + x2r;
  436. a[9] = x0i + x2i;
  437. a[12] = x2i - x0i;
  438. a[13] = x0r - x2r;
  439. x0r = x1r - x3i;
  440. x0i = x1i + x3r;
  441. a[10] = wk1r * (x0r - x0i);
  442. a[11] = wk1r * (x0r + x0i);
  443. x0r = x3i + x1r;
  444. x0i = x3r - x1i;
  445. a[14] = wk1r * (x0i - x0r);
  446. a[15] = wk1r * (x0i + x0r);
  447. k1 = 0;
  448. for (j = 16; j < n; j += 16) {
  449. k1 += 2;
  450. k2 = 2 * k1;
  451. wk2r = w[k1];
  452. wk2i = w[k1 + 1];
  453. wk1r = w[k2];
  454. wk1i = w[k2 + 1];
  455. wk3r = wk1r - 2 * wk2i * wk1i;
  456. wk3i = 2 * wk2i * wk1r - wk1i;
  457. x0r = a[j] + a[j + 2];
  458. x0i = a[j + 1] + a[j + 3];
  459. x1r = a[j] - a[j + 2];
  460. x1i = a[j + 1] - a[j + 3];
  461. x2r = a[j + 4] + a[j + 6];
  462. x2i = a[j + 5] + a[j + 7];
  463. x3r = a[j + 4] - a[j + 6];
  464. x3i = a[j + 5] - a[j + 7];
  465. a[j] = x0r + x2r;
  466. a[j + 1] = x0i + x2i;
  467. x0r -= x2r;
  468. x0i -= x2i;
  469. a[j + 4] = wk2r * x0r - wk2i * x0i;
  470. a[j + 5] = wk2r * x0i + wk2i * x0r;
  471. x0r = x1r - x3i;
  472. x0i = x1i + x3r;
  473. a[j + 2] = wk1r * x0r - wk1i * x0i;
  474. a[j + 3] = wk1r * x0i + wk1i * x0r;
  475. x0r = x1r + x3i;
  476. x0i = x1i - x3r;
  477. a[j + 6] = wk3r * x0r - wk3i * x0i;
  478. a[j + 7] = wk3r * x0i + wk3i * x0r;
  479. wk1r = w[k2 + 2];
  480. wk1i = w[k2 + 3];
  481. wk3r = wk1r - 2 * wk2r * wk1i;
  482. wk3i = 2 * wk2r * wk1r - wk1i;
  483. x0r = a[j + 8] + a[j + 10];
  484. x0i = a[j + 9] + a[j + 11];
  485. x1r = a[j + 8] - a[j + 10];
  486. x1i = a[j + 9] - a[j + 11];
  487. x2r = a[j + 12] + a[j + 14];
  488. x2i = a[j + 13] + a[j + 15];
  489. x3r = a[j + 12] - a[j + 14];
  490. x3i = a[j + 13] - a[j + 15];
  491. a[j + 8] = x0r + x2r;
  492. a[j + 9] = x0i + x2i;
  493. x0r -= x2r;
  494. x0i -= x2i;
  495. a[j + 12] = -wk2i * x0r - wk2r * x0i;
  496. a[j + 13] = -wk2i * x0i + wk2r * x0r;
  497. x0r = x1r - x3i;
  498. x0i = x1i + x3r;
  499. a[j + 10] = wk1r * x0r - wk1i * x0i;
  500. a[j + 11] = wk1r * x0i + wk1i * x0r;
  501. x0r = x1r + x3i;
  502. x0i = x1i - x3r;
  503. a[j + 14] = wk3r * x0r - wk3i * x0i;
  504. a[j + 15] = wk3r * x0i + wk3i * x0r;
  505. }
  506. }
  507. void cftmdl(int n, int l, fftType *a, fftType *w)
  508. {
  509. int j, j1, j2, j3, k, k1, k2, m, m2;
  510. fftType wk1r, wk1i, wk2r, wk2i, wk3r, wk3i;
  511. fftType x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i;
  512. m = l << 2;
  513. for (j = 0; j < l; j += 2) {
  514. j1 = j + l;
  515. j2 = j1 + l;
  516. j3 = j2 + l;
  517. x0r = a[j] + a[j1];
  518. x0i = a[j + 1] + a[j1 + 1];
  519. x1r = a[j] - a[j1];
  520. x1i = a[j + 1] - a[j1 + 1];
  521. x2r = a[j2] + a[j3];
  522. x2i = a[j2 + 1] + a[j3 + 1];
  523. x3r = a[j2] - a[j3];
  524. x3i = a[j2 + 1] - a[j3 + 1];
  525. a[j] = x0r + x2r;
  526. a[j + 1] = x0i + x2i;
  527. a[j2] = x0r - x2r;
  528. a[j2 + 1] = x0i - x2i;
  529. a[j1] = x1r - x3i;
  530. a[j1 + 1] = x1i + x3r;
  531. a[j3] = x1r + x3i;
  532. a[j3 + 1] = x1i - x3r;
  533. }
  534. wk1r = w[2];
  535. for (j = m; j < l + m; j += 2) {
  536. j1 = j + l;
  537. j2 = j1 + l;
  538. j3 = j2 + l;
  539. x0r = a[j] + a[j1];
  540. x0i = a[j + 1] + a[j1 + 1];
  541. x1r = a[j] - a[j1];
  542. x1i = a[j + 1] - a[j1 + 1];
  543. x2r = a[j2] + a[j3];
  544. x2i = a[j2 + 1] + a[j3 + 1];
  545. x3r = a[j2] - a[j3];
  546. x3i = a[j2 + 1] - a[j3 + 1];
  547. a[j] = x0r + x2r;
  548. a[j + 1] = x0i + x2i;
  549. a[j2] = x2i - x0i;
  550. a[j2 + 1] = x0r - x2r;
  551. x0r = x1r - x3i;
  552. x0i = x1i + x3r;
  553. a[j1] = wk1r * (x0r - x0i);
  554. a[j1 + 1] = wk1r * (x0r + x0i);
  555. x0r = x3i + x1r;
  556. x0i = x3r - x1i;
  557. a[j3] = wk1r * (x0i - x0r);
  558. a[j3 + 1] = wk1r * (x0i + x0r);
  559. }
  560. k1 = 0;
  561. m2 = 2 * m;
  562. for (k = m2; k < n; k += m2) {
  563. k1 += 2;
  564. k2 = 2 * k1;
  565. wk2r = w[k1];
  566. wk2i = w[k1 + 1];
  567. wk1r = w[k2];
  568. wk1i = w[k2 + 1];
  569. wk3r = wk1r - 2 * wk2i * wk1i;
  570. wk3i = 2 * wk2i * wk1r - wk1i;
  571. for (j = k; j < l + k; j += 2) {
  572. j1 = j + l;
  573. j2 = j1 + l;
  574. j3 = j2 + l;
  575. x0r = a[j] + a[j1];
  576. x0i = a[j + 1] + a[j1 + 1];
  577. x1r = a[j] - a[j1];
  578. x1i = a[j + 1] - a[j1 + 1];
  579. x2r = a[j2] + a[j3];
  580. x2i = a[j2 + 1] + a[j3 + 1];
  581. x3r = a[j2] - a[j3];
  582. x3i = a[j2 + 1] - a[j3 + 1];
  583. a[j] = x0r + x2r;
  584. a[j + 1] = x0i + x2i;
  585. x0r -= x2r;
  586. x0i -= x2i;
  587. a[j2] = wk2r * x0r - wk2i * x0i;
  588. a[j2 + 1] = wk2r * x0i + wk2i * x0r;
  589. x0r = x1r - x3i;
  590. x0i = x1i + x3r;
  591. a[j1] = wk1r * x0r - wk1i * x0i;
  592. a[j1 + 1] = wk1r * x0i + wk1i * x0r;
  593. x0r = x1r + x3i;
  594. x0i = x1i - x3r;
  595. a[j3] = wk3r * x0r - wk3i * x0i;
  596. a[j3 + 1] = wk3r * x0i + wk3i * x0r;
  597. }
  598. wk1r = w[k2 + 2];
  599. wk1i = w[k2 + 3];
  600. wk3r = wk1r - 2 * wk2r * wk1i;
  601. wk3i = 2 * wk2r * wk1r - wk1i;
  602. for (j = k + m; j < l + (k + m); j += 2) {
  603. j1 = j + l;
  604. j2 = j1 + l;
  605. j3 = j2 + l;
  606. x0r = a[j] + a[j1];
  607. x0i = a[j + 1] + a[j1 + 1];
  608. x1r = a[j] - a[j1];
  609. x1i = a[j + 1] - a[j1 + 1];
  610. x2r = a[j2] + a[j3];
  611. x2i = a[j2 + 1] + a[j3 + 1];
  612. x3r = a[j2] - a[j3];
  613. x3i = a[j2 + 1] - a[j3 + 1];
  614. a[j] = x0r + x2r;
  615. a[j + 1] = x0i + x2i;
  616. x0r -= x2r;
  617. x0i -= x2i;
  618. a[j2] = -wk2i * x0r - wk2r * x0i;
  619. a[j2 + 1] = -wk2i * x0i + wk2r * x0r;
  620. x0r = x1r - x3i;
  621. x0i = x1i + x3r;
  622. a[j1] = wk1r * x0r - wk1i * x0i;
  623. a[j1 + 1] = wk1r * x0i + wk1i * x0r;
  624. x0r = x1r + x3i;
  625. x0i = x1i - x3r;
  626. a[j3] = wk3r * x0r - wk3i * x0i;
  627. a[j3 + 1] = wk3r * x0i + wk3i * x0r;
  628. }
  629. }
  630. }
  631. void rftfsub(int n, fftType *a, int nc, fftType *c)
  632. {
  633. int j, k, kk, ks, m;
  634. fftType wkr, wki, xr, xi, yr, yi;
  635. m = n >> 1;
  636. ks = 2 * nc / m;
  637. kk = 0;
  638. for (j = 2; j < m; j += 2) {
  639. k = n - j;
  640. kk += ks;
  641. wkr = 0.5 - c[nc - kk];
  642. wki = c[kk];
  643. xr = a[j] - a[k];
  644. xi = a[j + 1] + a[k + 1];
  645. yr = wkr * xr - wki * xi;
  646. yi = wkr * xi + wki * xr;
  647. a[j] -= yr;
  648. a[j + 1] -= yi;
  649. a[k] += yr;
  650. a[k + 1] -= yi;
  651. }
  652. }
  653. void rftbsub(int n, fftType *a, int nc, fftType *c)
  654. {
  655. int j, k, kk, ks, m;
  656. fftType wkr, wki, xr, xi, yr, yi;
  657. a[1] = -a[1];
  658. m = n >> 1;
  659. ks = 2 * nc / m;
  660. kk = 0;
  661. for (j = 2; j < m; j += 2) {
  662. k = n - j;
  663. kk += ks;
  664. wkr = 0.5 - c[nc - kk];
  665. wki = c[kk];
  666. xr = a[j] - a[k];
  667. xi = a[j + 1] + a[k + 1];
  668. yr = wkr * xr + wki * xi;
  669. yi = wkr * xi - wki * xr;
  670. a[j] -= yr;
  671. a[j + 1] = yi - a[j + 1];
  672. a[k] += yr;
  673. a[k + 1] = yi - a[k + 1];
  674. }
  675. a[m + 1] = -a[m + 1];
  676. }

rdft4g.cpp : last modifiled at 2017/07/26 10:59:45(20,130 bytes)
created at 2016/04/11 11:17:20
The creation time of this html file is 2017/10/07 10:54:16 (Sat Oct 07 10:54:16 2017).